package mandel;

import java.awt.image.BufferedImage;
import java.awt.image.DataBufferInt;
import java.util.ArrayList;


public class Calculateur {
	
	// Taille découpage en X et en Y de la zone de dessin
	// Si trop petit : pas assez de granularité et donc de répartition sur l'ensemble des thread
	//  et moins de résultat rapide dans la zone centrale du zoom
	// Si trop grand : moins de possibilité d'avoir de grand rectangle fillés ou interpolés
	// L'optimum semble être environ 8
	private static final int SIZE_DECOUPAGE = 8; 
	private static int NB_DECOUPAGE = SIZE_DECOUPAGE*SIZE_DECOUPAGE;
	
	
	public static boolean bIscalculating;
	private static boolean bStopCalcul;
	public static int NB_ITER_MAX;
	
	private static Coordinates coords;
	private static OptimOption optims;

	private static double xToDrawFirst, yToDrawFirst;
	
	private static ColorRamp colorRamp;
	
	private static MainFrame mainFrame;
	private static ThreadPere threadPere;
	
	private static long debCalcul, finCalcul;
	
	private static int nbPixel_disp;
	private static int nbPixel_aa; // = TAILLE DES MATRIX (peut etre plus élevé si anti-aliasing)
	private static double[] matrix_val; // la ou vont etre stocké les resultat brut du calcul
	private static int[] matrix_color;// la ou vont etre stockées les couleurs
	private static int[] matrix_optim;// la ou vont etre stockées les info d'optim de calcul
	
	// qques stats
	private static int nbPixel_computed;
	private static int nbPixel_equationed; // déterminé comme étant dans le set via équation
	private static int nbPixel_filled;
	private static int nbPixel_interpolated;
	
	// Le PanelDessin fournit permet au Calculateur de forcer un repaint dès que le calcul est terminé
	public static void init(MainFrame mainFrame) {
		Calculateur.bIscalculating = false;
		Calculateur.bStopCalcul = false;

		Calculateur.mainFrame = mainFrame;
		
		Calculateur.coords = null;
		matrix_val = null;
		matrix_color = null;
	}
	

	// On est dans l'EDT
	// xToDrawFirst et yToDrawFirst sont les coordonnées de la souris dans le panel de dessin
	//  au moment de la demande de zoom : cela permet de commencer le calcul dans cette zone,
	//  qui est celle interessant l'utilisateur
	public static void startCalcul(Coordinates coords_new, int nbIterMax, double xToDrawFirst, double yToDrawFirst, ColorRamp colorRamp, OptimOption optims) {
		assert (bIscalculating == false);
		
		Calculateur.colorRamp = colorRamp;
		Calculateur.NB_ITER_MAX = nbIterMax;
		Calculateur.optims = optims;
		
		nbPixel_disp = coords_new.getNbPixels_disp();
		nbPixel_aa = coords_new.getNbPixels_aa();
				
		matrix_val = new double[nbPixel_aa];
		for (int i=0; i<nbPixel_aa; i++) matrix_val[i] = Double.NaN;
		
		// A priori, les bornes de calcul on changé
		// on initialise matrix_color avec le précédent calcul
		// ca permet à l'utilisateur d'avoir une premiere vision approchée du resultat final
		matrix_color = computePreImage_fromOld(matrix_color, Calculateur.coords, coords_new);
		
		matrix_optim = new int[nbPixel_aa];
		for (int i=0; i<nbPixel_aa; i++) matrix_optim[i] = ColorRamp.WHITE;
		
		// on conserve une copie, cela permet de faire un stretch de l'image en cas de zoom en guise d'attente
		// quand le coords géré par MainFrame est modifié avant un recalcul
		Calculateur.coords = coords_new.getCopy();
		
		// endroit ou il faut commencer le calcul
		Calculateur.xToDrawFirst = xToDrawFirst;
		Calculateur.yToDrawFirst = yToDrawFirst;
		
		mandel.Logger.log("---- DEBUT DU CALCUL ----");
		debCalcul = System.nanoTime();
		
		
		bIscalculating = true;
		bStopCalcul = false;
		
		threadPere = new ThreadPere();
		threadPere.start();
		
		// on renvoie des a present, pendant que le calcul tourne
	}
	
	
	// Ne retourne que quand le calcul est effectivement terminé
	public static void stopCalcul() {
		bStopCalcul = true;
		try {
			threadPere.join();
		} catch (InterruptedException ex) {}
		mandel.Logger.log("ARRET PREMATURE DU CALCUL !!");
	}
	
	
	// Appelé par le panel qui souhaite se redessiner, pendant le calcul ou bien à la fin
	// Remplit l'image avec l'ensemble des éléments déjà calculés
	public static void fillImageWithCurrentData(BufferedImage img, boolean bShowOptim) {
		// A-t-on deja lancé un calcul ?
		if (matrix_color == null)
			return;
		
		// Si la taille du panel a changé, on ne recopie rien
		if (img.getWidth() != coords.get_nbPixelLarg_disp() || img.getHeight() != coords.get_nbPixelHaut_disp())
			return;

		int[] imagePixelData = ((DataBufferInt)img.getRaster().getDataBuffer()).getData();
		
		int[] matrixToUse = bShowOptim ? matrix_optim : matrix_color;
		
		if (optims.tailleAA ==1) {
			// cas simple : recopie brute de la matrix_color
			System.arraycopy(matrixToUse, 0, imagePixelData, 0, nbPixel_aa);
		} else if (optims.tailleAA == 2) {
			// on fait une moyenne de 4 couleurs de matrix_color pour 1 pixel de l'image
			int nbPixelLarg_aa = coords.get_nbPixelLarg_aa();
			int nbPixelHaut_aa = coords.get_nbPixelHaut_aa();
			int indexDest = 0;
			int indexSrc = 0;

			for (int j=0; j<nbPixelHaut_aa; j+=2) {
				for (int i=0; i<nbPixelLarg_aa; i+=2) {
					indexSrc = j*nbPixelLarg_aa + i;
					int c1 = matrixToUse[indexSrc];
					indexSrc++;
					int c2 = matrixToUse[indexSrc];
					indexSrc+=nbPixelLarg_aa;
					int c3 = matrixToUse[indexSrc];
					indexSrc--;
					int c4 = matrixToUse[indexSrc];

					int newColor = ColorRamp.getMoyenneColor(c1, c2, c3, c4);
					imagePixelData[indexDest++] = newColor;
				}
			}
			
		} else if (optims.tailleAA == 3) {
			// on fait une moyenne de 9 couleurs de matrix_color pour 1 pixel de l'image
			int nbPixelLarg_aa = coords.get_nbPixelLarg_aa();
			int nbPixelHaut_aa = coords.get_nbPixelHaut_aa();
			int indexDest = 0;
			int indexSrc = 0;

			for (int j=0; j<nbPixelHaut_aa; j+=3) {
				for (int i=0; i<nbPixelLarg_aa; i+=3) {
					indexSrc = j*nbPixelLarg_aa + i;
					int c1 = matrixToUse[indexSrc];
					int c2 = matrixToUse[indexSrc+1];
					int c3 = matrixToUse[indexSrc+2];
					
					indexSrc = (j+1)*nbPixelLarg_aa + i;
					int c4 = matrixToUse[indexSrc];
					int c5 = matrixToUse[indexSrc+1];
					int c6 = matrixToUse[indexSrc+2];

					indexSrc = (j+2)*nbPixelLarg_aa + i;
					int c7 = matrixToUse[indexSrc];
					int c8 = matrixToUse[indexSrc+1];
					int c9 = matrixToUse[indexSrc+2];
					
					c1 = ColorRamp.getMoyenneColor(c1, c2, c3);
					c4 = ColorRamp.getMoyenneColor(c4, c5, c6);
					c7 = ColorRamp.getMoyenneColor(c7, c8, c9);
					
					int newColor = ColorRamp.getMoyenneColor(c1, c4, c7);
					imagePixelData[indexDest++] = newColor;
				}
			}

			// Plus lent !
//			for (int n=0; n<nbPixel_disp; n++) {
//				int c1 = matrix_color[indexSrc];
//				indexSrc++;
//				int c2 = matrix_color[indexSrc];
//				indexSrc+=nbPixelLarg_aa;
//				int c3 = matrix_color[indexSrc];
//				indexSrc--;
//				int c4 = matrix_color[indexSrc];
//				indexSrc = indexSrc - nbPixelLarg_aa + 2;
//				if (indexSrc%nbPixelLarg_aa == 0)
//					indexSrc+=nbPixelLarg_aa;
//				
//				int newColor = ColorRamp.getMoyenneColor(c1, c2, c3, c4);
//				imagePixelData[indexDest++] = newColor;
//			}
		}
		
		
		double fin = System.nanoTime();
	}
	
	
	private static int[] computePreImage_fromOld(int[] matrix_color_old, Coordinates coords_old, Coordinates coords_new) {
		int[] matrix_color_new = new int[nbPixel_aa];

		// Si on n'a jamais fait aucun calcul, alors on retourne une matrice vide
		if (matrix_color_old == null) {
			//for (int i=0; i<nbPixel_aa; i++) matrix_color_new[i] = ColorRamp.RED;
			return matrix_color_new;
		}
		// si changement de tailel de panel ou de taille d'AA, alors on retourne une matrice vide
		if (coords_old.get_nbPixelLarg_aa() != coords_new .get_nbPixelLarg_aa() || coords_old.get_nbPixelHaut_aa() != coords_new .get_nbPixelHaut_aa()) {
			//for (int i=0; i<nbPixel_aa; i++) matrix_color_new[i] = ColorRamp.RED;
			return matrix_color_new;
		}
		
		int nbPixelLarg_aa = coords_old.get_nbPixelLarg_aa();
		int nbPixelHaut_aa = coords_old.get_nbPixelHaut_aa();
		
		for (int j_new=0 ; j_new<nbPixelHaut_aa ; j_new++) {
			for (int i_new=0 ; i_new<nbPixelLarg_aa ; i_new++) {
				double x_new = coords_new.getXFromI_aa(i_new);
				double y_new = coords_new.getYFromJ_aa(j_new);
				int i_old = coords_old.getIFromX_aa(x_new);
				int j_old = coords_old.getJFromY_aa(y_new);
				
				int color;
				// il se peut que l'on tombe hors de l'ancienne matrice, si zoom arriere ou scroll
				if (i_old>=0 && i_old<nbPixelLarg_aa && j_old>=0 && j_old<nbPixelHaut_aa)
					color = matrix_color_old[j_old * nbPixelLarg_aa + i_old];
				else
					color = ColorRamp.WHITE;
				
				matrix_color_new[j_new*nbPixelLarg_aa + i_new] = color;
			}		
		}
		
		return matrix_color_new;
	}

	
	// Appelé par l'EDT
	// Les couleurs ont changé : il faut recalculer l'ensemble de matrix_color a partir de matrix_val
	static void changeColorRamp(ColorRamp colorRamp) {
		assert(bIscalculating == false);
		Calculateur.colorRamp = colorRamp;
		for (int i=0; i<nbPixel_aa; i++)
			matrix_color[i] = colorRamp.getColor(matrix_val[i]);
	}


	
	private static class ThreadPere extends Thread {
		private ArrayList<RectCalcul> vRectCalcul;
		private ArrayList<Integer> vIndexRectCalcul;
		private int indexNextRectToCompute;
		
		@Override
		public void run() {
			
			// On lance n thread sur différentes zones de calcul pour mieux répartir la charge de calcul
			//  et pour pouvoir commencer le calcul sur une endroit précis de l'image
			vRectCalcul = new ArrayList<RectCalcul>(NB_DECOUPAGE);
			vIndexRectCalcul = new ArrayList<Integer>(NB_DECOUPAGE);
			indexNextRectToCompute = 0;

			//******************************************************************
			// Création de différentes zones de calcul
			int largMatrix = coords.get_nbPixelLarg_aa();
			int hautMatrix = coords.get_nbPixelHaut_aa();
			for (int j=0; j<SIZE_DECOUPAGE; j++)
				for (int i=0; i<SIZE_DECOUPAGE; i++)
					vRectCalcul.add(new RectCalcul( (largMatrix*i)/SIZE_DECOUPAGE,
													(largMatrix*(i+1))/SIZE_DECOUPAGE-1,
													(hautMatrix*j)/SIZE_DECOUPAGE,
													(hautMatrix*(j+1))/SIZE_DECOUPAGE-1));
			
			//******************************************************************
			// Puis on va les lancer dans un certain ordre en partant du point (xToDrawFirst, yToDrawFirst)
			//  et en faisant des cercles concentriques
		
			// On commence par déterminer la case de départ
			int iPixel, jPixel;
			iPixel = coords.getIFromX_aa(xToDrawFirst);
			jPixel = coords.getJFromY_aa(yToDrawFirst);
			
			int iCase_centre, jCase_centre;
			iCase_centre = (int)( (double)iPixel / largMatrix * SIZE_DECOUPAGE );
			jCase_centre = (int)( (double)jPixel / hautMatrix * SIZE_DECOUPAGE );
			
			addRectCalcul(iCase_centre, jCase_centre);
			
			// Puis on fait des "cercles" concentriques sur cette case centrale
			//  jusqu'à parcourir toutes les cases
			for (int numCercle=1 ; numCercle<=SIZE_DECOUPAGE-1 ; numCercle++) {
				int tailleCercle = 2*numCercle + 1;
				
				// coin haut-gauche du "cercle"
				int iCase_hg = iCase_centre - numCercle;
				int jCase_hg = jCase_centre - numCercle;

				// bord haut
				for (int k=0; k<=tailleCercle-1; k++)
					addRectCalcul(iCase_hg+k, jCase_hg);
				// bord bas
				for (int k=0; k<=tailleCercle-1; k++)
					addRectCalcul(iCase_hg+k, jCase_hg+tailleCercle-1);
				// bord gauche
				for (int k=1; k<=tailleCercle-2; k++)
					addRectCalcul(iCase_hg, jCase_hg+k);
				// bord droit
				for (int k=1; k<=tailleCercle-2; k++)
					addRectCalcul(iCase_hg+tailleCercle-1, jCase_hg+k);

			}
			assert (vIndexRectCalcul.size()==NB_DECOUPAGE);
			
			//******************************************************************
			// Enfin, on peut créer et lancer les thread de calcul
			// Les n threads vont se répartir les rect de calcul via GiveMeJob();
			ThreadCalcul[] tabThread = new ThreadCalcul[optims.nbThread];
			for (int i=0; i<optims.nbThread; i++)
				tabThread[i] = new ThreadCalcul(this, i);

			for (int i=0; i<optims.nbThread; i++)
				tabThread[i].start();
			
			// On attend la fin de tous les threads de calcul
			try {
				for (int i=0; i<optims.nbThread; i++)
					tabThread[i].join();
			} catch (InterruptedException ex) {}
			
			// Calcul stats
			// Chaque thread doit faire ses propres incrementation de compteur perso
			//  aggrégé en fin de calcul. Si chaque thread incremente des compteur globaux
			//  ca buggue car incrémentation non atomique (il faudrait synchroniser l'incrementation ...)
			nbPixel_computed = 0;
			nbPixel_equationed = 0;
			nbPixel_filled = 0;
			nbPixel_interpolated = 0;
			for (int i=0; i<optims.nbThread; i++) {
				nbPixel_computed += tabThread[i].nbPixel_computed;
				nbPixel_equationed += tabThread[i].nbPixel_equationed;
				nbPixel_filled += tabThread[i].nbPixel_filled;
				nbPixel_interpolated += tabThread[i].nbPixel_interpolated;
			}
			
			bIscalculating = false;
			
			// On n'affiche la duree de calcul que si arret naturel apres l'ensemble du calcul
			if (!bStopCalcul) {
				finCalcul = System.nanoTime();
				Logger.log("---- FIN DU CALCUL : " + ( (double)(finCalcul-debCalcul) / 1000000.0) + " ms ----");
				if (nbPixel_aa != (nbPixel_computed+nbPixel_filled+nbPixel_interpolated))
					Logger.log("!!!!!!!!!!!!!!!! " + nbPixel_aa + " / " + (nbPixel_computed+nbPixel_filled+nbPixel_interpolated));
				Logger.log(" => Nb pixel calculés : " + UtilsCalcul.round( 100*((double)nbPixel_computed) / nbPixel_aa, 2) + " %");
				Logger.log("   -> Total                : " + nbPixel_aa);
				Logger.log("   -> Calculé complètement : " + nbPixel_computed);
				Logger.log("   -> Déterminé via equat° : " + nbPixel_equationed);
				Logger.log("   -> Rempliss. uniforme   : " + nbPixel_filled);
				Logger.log("   -> Rempliss. interpolé  : " + nbPixel_interpolated);
			}
			// Refresh du panel de dessin
			mainFrame.calculTermine();
		}
		
		// Fonction utile pour construire l'ordre des rect de calcul
		// i et j peuvent sortir du cadre, dans ce cas, on n'ajoute rien
		private void addRectCalcul (int iCase, int jCase) {
			if (iCase<0 || iCase >= SIZE_DECOUPAGE || jCase<0 || jCase >= SIZE_DECOUPAGE)
				return;
			vIndexRectCalcul.add(jCase * SIZE_DECOUPAGE + iCase);
		}
		
		// Appelé par les thread de calcul en soif de travail
		// Renvoie null quand tout est terminé
		synchronized protected RectCalcul giveMeJob() {
			if (indexNextRectToCompute == vRectCalcul.size())
				return null;
			int indexRect = vIndexRectCalcul.get(indexNextRectToCompute++);
			return vRectCalcul.get(indexRect);
		}
	}
	
	
	private static class ThreadCalcul extends Thread {
		ThreadPere papoun;
		int numThread;
		int nbPixelLarg_aa;
		
		// Stats propres à chaque thread, pour ensuite alimenter les stats globales
		int nbPixel_computed;
		int nbPixel_equationed;
		int nbPixel_filled;
		int nbPixel_interpolated;
		
		public ThreadCalcul(ThreadPere papoun, int numThread) {
			this.papoun = papoun;
			this.numThread = numThread;
			this.nbPixelLarg_aa = coords.nbPixelLarg_aa;
			this.nbPixel_computed = 0;
			this.nbPixel_filled = 0;
			this.nbPixel_interpolated = 0;
		}
			
		@Override
		public void run() {
			// Boucle sur les rect de calcul
			RectCalcul r;
			while(true) {
				r = papoun.giveMeJob();
				if (r==null)
					return; // plus rien à calculer !
				
				if (!optims.bOptimFill && !optims.bOptimInterpol)
					computeRect_withoutOptim(r.iMin, r.iMax, r.jMin, r.jMax);
				else
					computeRect_withOptim(r.iMin, r.iMax, r.jMin, r.jMax, true);
			}
		}
		
		
		// Calcul explicite de tous les pixel du rectangle sans les optims Fill et Interpol
		private void computeRect_withoutOptim(int iMin, int iMax, int jMin, int jMax) {
			double xC, yC;
			double val;
			// Boucle de calcul sur les pixels du rectangle de calcul
			for (int j=jMin ; j<=jMax ; j++) {
				yC = coords.getYFromJ_aa(j);
				for (int i=iMin ; i<=iMax ; i++) {
					// Arret prematuré du calcul par l'utilisateur
					if (bStopCalcul) 	return;

					xC = coords.getXFromI_aa(i);
					val = computePoint_double(xC, yC);

					matrix_val[j * nbPixelLarg_aa + i] = val;
					int color = colorRamp.getColor(val);
					matrix_color[j * nbPixelLarg_aa + i] = color;
//					matrix_optim[j * nbPixelLarg_aa + i] = color;
				}
			}
		}


		// Types d'optimisations possibles des rect
		private final int OPTIM_RECT_FILLABLE = 0;		// le top
		private final int OPTIM_RECT_INTERPOLABLE = 1;	// pas mal
		private final int OPTIM_RECT_NO = 2;			// bigre, on doit subdiviser le rectangle

		// Calcul avec l'une au moins des optims Fill ou Interpol telles que définies dans l'objet 'optims'
		private void computeRect_withOptim(int iMin, int iMax, int jMin, int jMax, boolean bDecoupVertical) {
			int i, j, offset;
			double xC, yC;
			double val;
			
			// Calcul des bords du rectangle
			// - bord haut
			j = jMin;
			yC = coords.getYFromJ_aa(j);
			for (i=iMin ; i<=iMax ; i++) {
				// Arret prematuré du calcul par l'utilisateur
				if (bStopCalcul) return;
				offset = j * nbPixelLarg_aa + i;
				// si deja calculé, on passe au pixel suivant
				if (!Double.isNaN(matrix_val[offset]))	continue;
				xC = coords.getXFromI_aa(i);
				val = computePoint_double(xC, yC);
				matrix_val[offset] = val;
				int color = colorRamp.getColor(val);
				matrix_color[offset] = color;
//				matrix_optim[offset] = color;
			}
			// - bord bas
			j = jMax;
			yC = coords.getYFromJ_aa(j);
			for (i=iMin ; i<=iMax ; i++) {
				// Arret prematuré du calcul par l'utilisateur
				if (bStopCalcul) return;
				offset = j * nbPixelLarg_aa + i;
				// si deja calculé, on passe au pixel suivant
				if (!Double.isNaN(matrix_val[offset]))	continue;
				xC = coords.getXFromI_aa(i);
				val = computePoint_double(xC, yC);
				matrix_val[offset] = val;
				int color = colorRamp.getColor(val);
				matrix_color[offset] = color;
//				matrix_optim[offset] = color;
			}
			// - bord gauche
			i = iMin;
			xC = coords.getXFromI_aa(i);
			for (j=jMin+1 ; j<=jMax-1 ; j++) {
				// Arret prematuré du calcul par l'utilisateur
				if (bStopCalcul) return;
				offset = j * nbPixelLarg_aa + i;
				// si deja calculé, on passe au pixel suivant
				if (!Double.isNaN(matrix_val[offset]))	continue;
				yC = coords.getYFromJ_aa(j);
				val = computePoint_double(xC, yC);
				matrix_val[offset] = val;
				int color = colorRamp.getColor(val);
				matrix_color[offset] = color;
//				matrix_optim[offset] = color;
			}
			// - bord droit
			i = iMax;
			xC = coords.getXFromI_aa(i);
			for (j=jMin+1 ; j<=jMax-1 ; j++) {
				// Arret prematuré du calcul par l'utilisateur
				if (bStopCalcul) return;
				offset = j * nbPixelLarg_aa + i;
				// si deja calculé, on passe au pixel suivant
				if (!Double.isNaN(matrix_val[offset]))	continue;
				yC = coords.getYFromJ_aa(j);
				val = computePoint_double(xC, yC);
				matrix_val[offset] = val;
				int color = colorRamp.getColor(val);
				matrix_color[offset] = color;
//				matrix_optim[offset] = color;
			}
			
			// Si le rectangle n'est constitué que d'un bord, on a déjà tout parcouru !
			if (iMax-iMin+1 <=2 || jMax-jMin+1 <=2) {
				return;
			}
			
			// Sinon, on teste quel type d'optimisation est possible
			int typeOptim = getPossibleOptimisation(iMin, iMax, jMin, jMax);
			if (typeOptim == OPTIM_RECT_FILLABLE)
				fillInteriorRect(iMin, iMax, jMin, jMax);
			else if (typeOptim == OPTIM_RECT_INTERPOLABLE)
				interpoleInteriorRect(iMin, iMax, jMin, jMax);
			else {
				// Sinon, pas d'optimisation possible : on découpe le rectangle en 2 puis appel recursif
				// Alternativement, on découpe le rectangle verticalement, puis horizontalement
				if (bDecoupVertical) {
					computeRect_withOptim(iMin, (iMax+iMin)/2, jMin, jMax, !bDecoupVertical);
					computeRect_withOptim ((iMax+iMin)/2, iMax, jMin, jMax, !bDecoupVertical);					
				} else {
					computeRect_withOptim (iMin, iMax, jMin, (jMax+jMin)/2, !bDecoupVertical);
					computeRect_withOptim (iMin, iMax, (jMax+jMin)/2, jMax, !bDecoupVertical);
				}

			}
		}

		
		private int getPossibleOptimisation (int iMin, int iMax, int jMin, int jMax) {
			boolean bIsSameColor = true;
			int i, j;
			int deltaColor;
			
			// - bord haut
			for (i=iMin ; i<=iMax-1 ; i++) {
				deltaColor = getDeltaColor (jMin * nbPixelLarg_aa + i, jMin * nbPixelLarg_aa + i+1);
				if (deltaColor>optims.seuilMaxInterpol) return OPTIM_RECT_NO; // aucune optim possible
				if (deltaColor>optims.seuilMaxFill) bIsSameColor = false; // OPTIM_RECT_FILLABLE nous est désormais interdit
			}
			// - bord bas
			for (i=iMin ; i<=iMax-1 ; i++) {
				deltaColor = getDeltaColor (jMax * nbPixelLarg_aa + i, jMax * nbPixelLarg_aa + i+1);
				if (deltaColor>optims.seuilMaxInterpol) return OPTIM_RECT_NO; // aucune optim possible
				if (deltaColor>optims.seuilMaxFill) bIsSameColor = false; // OPTIM_RECT_FILLABLE nous est désormais interdit
			}
			// - bord gauche
			for (j=jMin ; j<=jMax-1 ; j++) {
				deltaColor = getDeltaColor (j * nbPixelLarg_aa + iMin, (j+1) * nbPixelLarg_aa + iMin);
				if (deltaColor>optims.seuilMaxInterpol) return OPTIM_RECT_NO; // aucune optim possible
				if (deltaColor>optims.seuilMaxFill) bIsSameColor = false; // OPTIM_RECT_FILLABLE nous est désormais interdit
			}
			// - bord droit
			for (j=jMin ; j<=jMax-1 ; j++) {
				deltaColor = getDeltaColor (j * nbPixelLarg_aa + iMax, (j+1) * nbPixelLarg_aa + iMax);
				if (deltaColor>optims.seuilMaxInterpol) return OPTIM_RECT_NO; // aucune optim possible
				if (deltaColor>optims.seuilMaxFill) bIsSameColor = false; // OPTIM_RECT_FILLABLE nous est désormais interdit
			}
			// si on arrive la, c'est que l'on est soit INTERPOLABLE soit FILLABLE
			if (bIsSameColor)
				return OPTIM_RECT_FILLABLE;
			else
				return OPTIM_RECT_INTERPOLABLE;
		}
		
		
		// A partir de 2 pixels adjacents d'un bord de Rect, détermine l'écart de couleur
		private int getDeltaColor(int offset1, int offset2) {
			int color1 = matrix_color[offset1];
			int color2 = matrix_color[offset2];
			return ColorRamp.getDeltaColor(color1, color2);
		}
		

		// La couleur de remplissage est celle du bord
		private void fillInteriorRect(int iMin, int iMax, int jMin, int jMax) {
			if (!optims.bOptimFill)
				// On aurait droit au Fill mais l'utilisateur l'a désactivé
				computeRect_withoutOptim(iMin+1, iMax-1, jMin+1, jMax-1);
			else {
				int offset;
				// coin haut gauche
				offset = jMin * nbPixelLarg_aa + iMin;
				double val_hg = matrix_val[offset];
				int color_hg = matrix_color[offset];
				for (int j=jMin+1 ; j<=jMax-1 ; j++) {
					for (int i=iMin+1 ; i<=iMax-1 ; i++) {
						this.nbPixel_filled++;
						offset = j * nbPixelLarg_aa + i;
						matrix_val[offset] = val_hg;
						matrix_color[offset] = color_hg;
						matrix_optim[offset] = ColorRamp.GREEN_LIGHT;
						
					}
				}
			}
		}
		
		private void interpoleInteriorRect(int iMin, int iMax, int jMin, int jMax) {
			if (!optims.bOptimInterpol)
				// On aurait droit au Fill mais l'utilisateur l'a désactivé
				computeRect_withoutOptim(iMin+1, iMax-1, jMin+1, jMax-1);
			else {
				int offset;
				// coin haut-gauche
//				offset = jMin * nbPixelLarg_aa + iMin;
//				double iter_hg = matrix_iter[offset];
//				int color_hg = matrix_color[offset];
				for (int j=jMin+1 ; j<=jMax-1 ; j++) {
					for (int i=iMin+1 ; i<=iMax-1 ; i++) {
						this.nbPixel_interpolated++;
						
						// En fait, on interpole les nbIter et on calcule la couleur ensuite
						// Si on interpole les couleurs, on ne saura pas quelle iter mettre dans matrix_iter !
						double valHaut = matrix_val [jMin * nbPixelLarg_aa + i];
						double valBas = matrix_val [jMax * nbPixelLarg_aa + i];
						double interpoleVertical = (valHaut*(jMax-j) + valBas*(j-jMin)) / (jMax-jMin);
						
						double valGauche = matrix_val [j * nbPixelLarg_aa + iMin];
						double valDroite = matrix_val [j * nbPixelLarg_aa + iMax];
						double interpoleHorizontal = (valGauche*(iMax-i) + valDroite*(i-iMin)) / (iMax-iMin);
						
						double valInterpol = (interpoleVertical + interpoleHorizontal)/2;
						
						offset = j * nbPixelLarg_aa + i;
						matrix_val [offset] = valInterpol;
						matrix_color[offset] = colorRamp.getColor(valInterpol);
						matrix_optim[offset] = ColorRamp.RED_LIGHT;
					}
				}
			}
		}
		
		
		// Le coeur du calcul !! La ou se passe 99% du temps de calcul
		private double computePoint_double(double xC, double yC) {
			return computePoint_double_DE(xC, yC);
		}
			
		// Renvoie l'escape time (=nb d'itération pour dépasser le BailOut)
		private double computePoint_double_iter(double xC, double yC) {
			
			// Test simple d'appartenance à l'ensemble
			double q = (xC-0.25) * (xC-0.25) + yC*yC;
			if ( (xC+1)*(xC+1) + yC*yC < 0.0625 ||
					q*(q+(xC-0.25)) < yC*yC/4  ) {
				// on est sur d'etre dans l'ensemble
				this.nbPixel_equationed++;
				return Double.POSITIVE_INFINITY;

			} else {
				// Calcul complet
				this.nbPixel_computed++;

				double xZ = xC;
				double yZ = yC;
				double xZ2;
				double yZ2;

				int numIter = 0;
				while(numIter<NB_ITER_MAX) {
					xZ2 = xZ*xZ;
					yZ2 = yZ*yZ;
					yZ = 2*xZ*yZ + yC;
					xZ = xZ2 - yZ2 + xC;
					if (xZ2 + yZ2 > 4) return numIter;
									
					numIter++;
				}
				return Double.POSITIVE_INFINITY;
			}
		}
		
		// Renvoie la distance estimate (-1 signifie dans l'ensemble)
		private double computePoint_double_DE(double xC, double yC) {
			
			// Test simple d'appartenance à l'ensemble
			double q = (xC-0.25) * (xC-0.25) + yC*yC;
			if ( (xC+1)*(xC+1) + yC*yC < 0.0625 ||
					q*(q+(xC-0.25)) < yC*yC/4  ) {
				// on est sur d'etre dans l'ensemble
				this.nbPixel_equationed++;
				return -1;

			} else {
				// Calcul complet
				this.nbPixel_computed++;

				double xZ = xC;
				double yZ = yC;
				// données temp pour optim
				double xZ2;
				double yZ2;
				// dérivées mise à jour à chaque itération calculer la distance estimate en fin de calcul
				double xDZ=1, xDZ_temp;
				double yDZ=0;

				int numIter = 0;
				while(numIter<NB_ITER_MAX) {
					
					// Calcul Z
					xZ2 = xZ*xZ;
					yZ2 = yZ*yZ;
					yZ = 2*xZ*yZ + yC;
					xZ = xZ2 - yZ2 + xC;
					
					// Calcul Z' (Z' = 2.Z.Z' + 1)
					xDZ_temp = 2 * (xZ*xDZ - yZ*yDZ) + 1;
					yDZ      = 2 * (yZ*xDZ + xZ*yDZ);
					xDZ = xDZ_temp;
					
					if (xZ2 + yZ2 > 4) {
						// calcul du distance estimate
						// DE = |Z| * ln(|Z|) / |Z'|
						double normZ  = Math.sqrt(xZ *xZ   + yZ *yZ);
						double normDZ = Math.sqrt(xDZ*xDZ  + yDZ*yDZ);
						return normZ * Math.log(normZ) / normDZ;
					}
									
					numIter++;
				}
				return -1;
			}
		}
	}
	
	
	
	private static class RectCalcul {
		int iMin, iMax, jMin, jMax;

		public RectCalcul(int iMin, int iMax, int jMin, int jMax) {
			this.iMin = iMin;
			this.iMax = iMax;
			this.jMin = jMin;
			this.jMax = jMax;
		}
	}

	
	/* Données empiriques, avec DELTA_X_AT_ZOOM_1 = 4, largeur ecran=600
		zoom		taille pixel	itermin
		1			0,006666667		200
		30			0,000222222		500
		3900		1,7094E-06		2000
		65000		1,02564E-07		5000
		200000		3,33333E-08		8000
		4,00E+06	1,66667E-09		50000
		2,00E+07	3,33333E-10		100000
		5,00E+08	1,33333E-11		350000

		Ce qui donne : ln(itermin) = 0.373 * ln(1/taillePixel) + 5.3
	 */
	static int getNbIterFromZoom(double zoomLevel, int nbPixelLarg_aa) {
		double taillePixel = (Coordinates.DELTA_X_AT_ZOOM_1/zoomLevel)/nbPixelLarg_aa;
		return (int) (Math.exp( 0.373 * Math.log(1/taillePixel) + 3.4 ));
	}
}
